//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   CartMesh.hh
 * \author Jeremy Roberts
 * \date   07/06/2011
 * \brief  Mesh member definitions.
 * \note   Copyright (C) Jeremy Roberts.
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef CARTMESH_HH_
#define CARTMESH_HH_

#include <vector>
#include <map>

#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "../utilities/DBC.hh"
#include "../utilities/SP.hh"
#include "../utilities/Soft_Equivalence.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class Base_Mesh
 * \brief Abstract base mesh class.
 *
 * We use an abstract base mesh that helps avoid using the templated Mesh as
 * a method argument.
 *
 * \todo I think many of the arguments for Base_Mesh and Mesh methods can
 *       and should be "const" following Meyers.
 *
 * \todo Rename to BaseMesh.
 *
 * \example test/testMesh.cc
 *
 */
//===========================================================================//
class Base_Mesh
{

public:

  //@{
  //! Useful typedefs for classes that use mesh
  typedef typename util::SP<Base_Mesh>  SP_mesh;
  typedef std::vector<double>           Vec_Dbl;
  typedef std::vector<int>              Vec_Int;
  typedef std::vector<int>              Vec_Id;
  //@}

  Base_Mesh(int number_cells, int number_faces = 0) :
    d_number_cells(number_cells),
    d_number_faces(number_faces)
  {}

  virtual ~Base_Mesh() = 0;

  // Return the number of cells
  int number_cells() const
  {
    return d_number_cells;
  }

  // Return the number of faces
  int number_faces() const
  {
    return d_number_faces;
  }

  // Return the number of cells in a given direction
  virtual int number_cells_direction(int dimension) const = 0;

  // Convert (i,j,k) to a cardinal index.
  virtual int convert(int i, int j = 0, int k = 0) const = 0;

  // Return fine cell width along \e (i,j,k).
  virtual double width(int ijk, int dimension = 0) const = 0;

  // Return inverse of fine cell width along \e (i,j,k).
  virtual double inv_width(int ijk, int dimension = 0) const = 0;

  // Cell center along \e (i,j,k)
  virtual double center(int ijk, int dimension = 0) const = 0;

  // Cell edges along a given direction.
  virtual const Vec_Dbl& edges(int dimension = 0) const = 0;

  // Region id.
  virtual int region_index(int i) const = 0;

  // Set cell id.
  virtual void set_cell_id(const int key, const Vec_Id &value) = 0;

  // Get cell id.
  virtual const Vec_Id& get_cell_id(const int key) = 0;

  // Nicely print me
  virtual void display() = 0;

protected:

  // Total number of cells in this mesh.
  const int d_number_cells;

  // Total number of faces in this mesh.
  const int d_number_faces;

};

//===========================================================================//
/*!
 * \class Mesh
 * \brief 1D/2D/3D orthogonal mesh for SLABTRAN
 *
 * This is a templated class for building a 1D/2D/3D Cartesian mesh.
 *
 * The mesh can be defined in one of three ways:
 * - defined coarse mesh edges in which a set number of fine meshes is defined
 * - defined fine mesh edges
 * - a single total mesh width (from origin) and mesh count for all dimensions
 *
 * Each mesh has a unique "region" id.  For meshes defined by coarse
 * meshes, each coarse region is given an id, and each and each fine mesh
 * inherits the appropriate id.  For a mesh defined by fine meshes directly
 * the cell id's are just their cardinal index.  For the single mesh width
 * option, their is just one region id.
 *
 * Using a region id facilitates material placement, among other things.
 * All ordering is defined along \f$x\f$, then \f$y\f$, then \f$z\f$.
 *
 * Note, use is made of "dimensioned vectors" based on the
 * \ref DimVecTraits class.  It helps facilitate generic programming.  For
 * example, suppose we wish to store cell edges in a generic container
 * of vectors, the type of which is defined as "Cell_Edges".  For 1d,
 * this vector should be just that, a vector<double>.  For 2d and 3d,
 * we want a vector of these, vector<double>, for which we use
 * \ref Vector_Lite.  Because these typedef's depend on the template
 * parameter, we use the DimVecTraits to define for us the appropriate
 * typedef.  2d and 3d are handled the same way, while 1d get's a
 * special treatment, as there is no use having a Vector_Lite of size
 * one in which something else resides.
 *
 * \todo can we use c++ template metaprogramming to eliminate repetitious
 *       loops?
 * \todo test each constructor and public accessor with established
 *       hardcoded examples
 * \todo verify monotonic increasing edges -- checked when returning width
 * \todo add a derived class for generating cmfd grids
 * \todo some way to designate surfaces for keeping angular fluxes
 * \todo implement method for face counting (not as trivial as cell counts)
 * \todo It would be highly useful to develop a modeling scheme in which
 *       a user defines pin cells, assemblies of pins, and cores of
 *       assemblies.  For reactor problems, or main purpose, we will need
 *       pin powers.  Having pins and assemblies be uniquely identified would
 *       be very helpful.  This shouldn't change the underlying discretization
 *       as implemented in \ref Mesh.  However, we need a sensible way to
 *       store "pin", "assembly", "region", "material", and any other
 *       identifiers that may be associated with each fine mesh.
 * \todo after a bit of thought, i think mesh should only be constructed
 *       via the fine mesh edges.  Going from temporaries (e.g. coarse
 *       mesh edges) to fine should be handled by an input processor.  The
 *       meshid's could then be assigned at the same time (after having
 *       processed the geometry description).  Hence, I should 1. remove
 *       all but the needed constructor, and 2. remove the Region index
 *       routines.
 */
 /*!
 * \example test/testMesh.cc
 *
 * Test of Mesh.
 */
//===========================================================================//

template <class D>
class Mesh : public Base_Mesh
{

public:

  /// \name Typedefs
  //\{
  /// Smart pointer to me
  typedef typename util::SP<Mesh<D> > SP_mesh;
  /// My base
  typedef Base_Mesh Base;
  /// A double vector
  typedef Base::Vec_Dbl Vec_Dbl;
  /// An integer vector
  typedef Base::Vec_Int Vec_Int;
  /// A vector of id's (just an integer vector)
  typedef Base::Vec_Id Vec_Id;
  /// A "dimensioned" integer (for one count in each direction)
  typedef typename DimVectorTraits<int, D>::DimVector Dim_Vector;
  /// A "dimensioned" integer vector (for several counts in each direction)
  typedef typename DimVectorTraits<Vec_Int, D>::DimVector Count_Vector;
  /// A "dimensioned" double (for one value in each direction)
  typedef typename DimVectorTraits<double, D>::DimVector Space_Vector;
  /// A "dimensioned" double vector (for several values in each direction)
  typedef typename DimVectorTraits<Vec_Dbl, D>::DimVector Cell_Edges;
  /// A "dimensioned" double vector (for several values in each direction)
  typedef typename DimVectorTraits<Vec_Dbl, D>::DimVector Cell_Widths;
  /// Cell id map
  typedef typename std::map<const int, Vec_Id>  Cell_Ids;
  /// Cell id map value type
  typedef Cell_Ids::value_type                  Cell_Ids_Value_Type;
  //\}

  // Public Members

  /*
   * \brief Constructor using coarse edges and fine counts
   *
   * \param coarse_region_edges   edges of coarse meshes
   * \param fine_per_coarse       number fine meshes in each coarse
   */
  explicit Mesh(Cell_Edges &coarse_region_edges,
                Count_Vector &fine_per_coarse);

  /*
   * \brief Constructor using (fine) mesh edges
   *
   * \param edges   edges of meshes
   */
  explicit Mesh(Cell_Edges &edges);

  /*
   * \brief Constructor using one width and fine count for all directions
   *
   * \param mesh_width   width of mesh mesh in each direction
   * \param fine_count   number of fine meshes in each direction
   */
  explicit Mesh(double mesh_width, int fine_count);

  // Default destructor
  ~Mesh(){};

  // >>> PUBLIC INTERFACE

  // INLINE

  //! \name Accessors
  //\{
  // Number of cells in a direction
  int number_cells_direction(int dimension = 0) const;

  // Convert (i,j,k) to a cardinal index.
  int convert(int i, int j = 0, int k = 0) const;

  // Return fine cell width along \e (i,j,k).
  double width(int ijk, int dimension = 0) const;

  // Return inverse of fine cell width along \e (i,j,k).
  double inv_width(int ijk, int dimension = 0) const;

  // Cell center along \e (i,j,k)
  double center(int ijk, int dimension = 0) const;

  // Cell edges along a given direction.
  const Vec_Dbl& edges(int dimension = 0) const
  {
    return d_edges[dimension];
  }

  /// Region id.
  int region_index(int i) const;

  // Set cell id.

  /*
   * \brief   Set cell id.
   * \param   key    name of the id (enum type)
   * \param   value  vector of cell id's
   */
  void set_cell_id(const int key, const Vec_Id &value);

  /*
   * \brief   Get cell id.
   *
   * This should \em not let the client receive nothing, i.e
   * we should bomb if we don't have it. (Or, we should
   * come up with a situation where we'd want conditional
   * treatment.)
   *
   * Also, since we incur a small cost for "finding" the
   * right data, it is highly recommended that the client
   * establish a reference once and use it throughout if
   * possible.  Local indexing can then be done via the
   * usual \ref Mesh::convert approach.
   *
   * \param   key    name of the id (enum type)
   * \return         vector of cell id's
   */
  const Vec_Id& get_cell_id(const int key);

  //\}

  /// Display the mesh edges and other information.
  void display();

protected:

  //---------------------------------------------------------------------------//
  // \name CONST DATA
  //\{
  /// Spatial dimension
  const int d_dimension;
  /// Number of cells in each dimension in this mesh.
  const Dim_Vector d_number_cells_direction;
  //\}

  //---------------------------------------------------------------------------//
  // \name MEMBER DATA
  //\{

  /// Fine mesh edges.
  Cell_Edges d_edges;

  /// Fine mesh widths.
  Cell_Widths d_width;

  /// Inverse fine mesh widths.
  Cell_Widths d_inv_width;

  /// Region id's
  Vec_Int d_region_index;

  /*!
   *  \brief Cell id's
   *
   *  One element of the cell id map contains a key (some integer from
   *  \ref Cell_Id_Types) and an integer vector that assigns one value per
   *  mesh cell.  These cell id's are used for material assignments,
   *  region associations (e.g. pin, assembly, fuel/mod.), etc.  Strings
   *  could alternatively be used as keys, but it's easier to add a
   *  needed id type to an enum list (so that trying to grab an invalid
   *  id type is a compile time error).  Using this structure allows
   *  us to add only the id's specified by the user (rather than
   *  having to modify the class for any additions).
   *
   *  For reference, std::map is used as follows:
   *  \code
   *    // generate e.g. Vec_Int matid and insert with key IDTYPE
   *    d_cell_ids.insert(d_cell_ids::value_type(MATERIAL_ID, matid));
   *    // retrieve matid
   *    const Vec_Id &matidref = d_cell_ids.find(MATERIAL_ID)->second;
   *  \endcode
   *
   */
  Cell_Ids d_cell_ids;

  //\}

  //---------------------------------------------------------------------------//
  // >>> Implementation

  /*
   * \brief   Compute total number of cells.
   * \param   V counts per coarse mesh in each direction
   * \return  total number of cells
   */
  int CountAllCells(Count_Vector &V);

  /*
   * \brief   Compute total number of cells (alternative)
   * \param   V fine edges in each dimension
   * \return  total number of cells
   */
  int CountAllCells(Cell_Edges &V);

  /*
   * \brief   Compute total number of cells in each dimension
   * \param   V counts per coarse mesh in each direction
   * \return  vector of counts in each column
   */
  Dim_Vector CountDimCells(Count_Vector &V);

  /*
   * \brief   Compute total number of cells in each dimension
   * \param   V fine mesh edges
   * \return  vector of counts in each column
   */
  Dim_Vector CountDimCells(typename Mesh<D>::Cell_Edges &V);

  /*
   * \brief   Generate region index from coarse mesh counts
   * \param   V counts per coarse mesh in each direction
   */
  void GenerateRegionIndex(Count_Vector &V);
};

//---------------------------------------------------------------------------//
// SPECIALIZATION DECLARATIONS
//---------------------------------------------------------------------------//

// Specialization for 1d.
template <>
int Mesh<_1D>::number_cells_direction(int dimension) const;

// Return edges: specialization to 1d
template <>
const Mesh<_1D>::Vec_Dbl& Mesh<_1D>::edges(int d) const;

// Specialization for 1d.
template <>
int Mesh<_1D>::CountAllCells(Mesh<_1D>::Count_Vector &V);

// Specialization for 1d.
template <>
Mesh<_1D>::Dim_Vector
Mesh<_1D>::CountDimCells(Mesh<_1D>::Count_Vector &V);

// Specialization for 1d.
template <>
int Mesh<_1D>::CountAllCells(Mesh<_1D>::Cell_Edges &V);

// Specialization for 1d.
template <>
Mesh<_1D>::Dim_Vector
Mesh<_1D>::CountDimCells(Mesh<_1D>::Cell_Edges &V);

template <>
void Mesh<_2D>::GenerateRegionIndex
(typename Mesh<_2D>::Count_Vector &fine_per_coarse);

template <>
void Mesh<_1D>::GenerateRegionIndex
(typename Mesh<_1D>::Count_Vector &fine_per_coarse);

} // end namespace slabtran

//---------------------------------------------------------------------------//
// INLINE FUNCTIONS
//---------------------------------------------------------------------------//

#include "tran/CartMesh.i.hh"

#endif // CARTMESH_HH_

//---------------------------------------------------------------------------//
//              end of Mesh.hh
//---------------------------------------------------------------------------//
